
. /* ------------------------------------------------------------------------ ** 
>     C. Wunder, A. Wiencierz, J. Schwarze, H. Küchenhoff:    
>     Well-being over the life span: semiparametric evidence from British 
>     and German longitudinal data
>     
>     ESTIMATION OF SEMIPARAMETRIC REGRESSIONS
> 
>     Data source:        British Household Panel Survey (BHPS), years 1996-2008.
>     Data organization:  person-year observations
> ** ------------------------------------------------------------------------ */ 
. 
. /*  --------( load data )-------------------------------------------------- */
. version 10

. matrix  drop _all

. use     BHPS.dta, clear
(PanelWhiz-BHPS: [SWB_and_Age_1991-2008]       [23 Jul 2009/17:00:36]  john@panel)

. rename age x

. rename gls_ls gls_y

. 
. /*  --------( set knots )-------------------------------------------------- */
. mkspline    spline 16 = x, marginal displayknots pctile

             |     knot1      knot2      knot3      knot4      knot5      knot6      knot7 
-------------+-----------------------------------------------------------------------------
           x |        21         25         29         32         35         38         41 

             |     knot8      knot9     knot10     knot11     knot12     knot13     knot14 
-------------+-----------------------------------------------------------------------------
           x |        44         48         52         56         60         65         70 

             |    knot15 
-------------+-----------
           x |        77 

. global      N_knots = r(N_knots)            

. global      N_spl = r(N_knots)+1            

. matrix      knots = r(knots)'

. 
. /*  --------( define macros & prepare data )------------------------------ */
. global exo_vars ///
>          ln_hhinc       ln_hhnpers      female          bad_health      /// 
>          educ_mid       educ_high       jstat_ausb      jstat_empl      ///
>          jstat_unem     jstat_reti      mstat2          mstat3          ///
>          mstat4         mstat5          mstat6          attr_in_1       ///
>          attr_in_2      attr_in_3       period2         period3         ///
>          period4        period5         period6         period7         /// 
>          period8        period9         period10        one

. 
. tokenize ${exo_vars}

. 
. global  gls_exo_vars ///
>         gls_`1'  gls_`2'  gls_`3'  gls_`4'  gls_`5'  gls_`6'  gls_`7'  gls_`8'  ///
>         gls_`9'  gls_`10' gls_`11' gls_`12' gls_`13' gls_`14' gls_`15' gls_`16' ///
>         gls_`17' gls_`18' gls_`19' gls_`20' gls_`21' gls_`22' gls_`23' gls_`24' ///
>         gls_`25' gls_`26' gls_`27' gls_`28'

. 
. global copy_gls_exo_vars ///
>         _gls_`1'  _gls_`2'  _gls_`3'  _gls_`4'  _gls_`5'  _gls_`6'  _gls_`7'  _gls_`8'  /
> //
>         _gls_`9'  _gls_`10' _gls_`11' _gls_`12' _gls_`13' _gls_`14' _gls_`15' _gls_`16' /
> //
>         _gls_`17' _gls_`18' _gls_`19' _gls_`20' _gls_`21' _gls_`22' _gls_`23' _gls_`24' /
> //
>         _gls_`25' _gls_`26' _gls_`27' _gls_`28' 

. 
. global TPF ///  
>         spline2     spline3     spline4     spline5     spline6     spline7     ///
>         spline8     spline9     spline10    spline11    spline12    spline13    ///
>         spline14    spline15    spline16  

. 
. foreach var of global TPF {
  2.     qui replace `var' = `var'^3
  3. }

. qui gen spline1_sq = spline1^2

. qui gen spline1_cu = spline1^3

. 
. global spline_vars spline1 spline1_sq spline1_cu $TPF

. foreach var of global spline_vars{
  2.     qui bysort pid: egen im_`var' = mean(`var')
  3.     qui gen gls_`var' = `var' - theta_i*im_`var'
  4. }

. drop im_*

. foreach var of global spline_vars{
  2.     qui gen _gls_`var' = gls_`var'
  3. }

. 
. /*  --------( estimation of semiparametric model eq. in 4 )----------------- */
. xtmixed gls_y gls_spline1 gls_spline1_sq gls_spline1_cu $gls_exo_vars, noconst /// 
>         || _all: gls_spline2-gls_spline$N_spl, noconstant cov(identity) ///
>         , technique(bfgs) emdots

Performing EM optimization: ....................

Performing gradient-based optimization: 

Iteration 0:   log restricted-likelihood = -162040.97  
Iteration 1:   log restricted-likelihood = -162040.97  
Iteration 2:   log restricted-likelihood = -162040.97  

Computing standard errors:

Mixed-effects REML regression                   Number of obs      =    123656
Group variable: _all                            Number of groups   =         1

                                                Obs per group: min =    123656
                                                               avg =  123656.0
                                                               max =    123656


                                                Wald chi2(31)      = 196437.05
Log restricted-likelihood = -162040.97          Prob > chi2        =    0.0000

--------------------------------------------------------------------------------
         gls_y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
   gls_spline1 |  -.3427377   .1959375    -1.75   0.080    -.7267682    .0412929
gls_spline1_sq |   .0114871   .0086841     1.32   0.186    -.0055334    .0285077
gls_spline1_cu |  -.0001306   .0001294    -1.01   0.313    -.0003841     .000123
  gls_ln_hhinc |   .0765989   .0065047    11.78   0.000     .0638499    .0893479
gls_ln_hhnpers |  -.0861881   .0119343    -7.22   0.000    -.1095789   -.0627973
    gls_female |   .0703592   .0133224     5.28   0.000     .0442479    .0964706
gls_bad_health |  -.2450499   .0080089   -30.60   0.000    -.2607471   -.2293528
  gls_educ_mid |   .0156623     .01705     0.92   0.358    -.0177552    .0490797
 gls_educ_high |    .033674   .0174098     1.93   0.053    -.0004485    .0677966
gls_jstat_ausb |    .328279   .0236117    13.90   0.000      .282001     .374557
gls_jstat_empl |   .2414269   .0122063    19.78   0.000      .217503    .2653507
gls_jstat_unem |   -.129368   .0199301    -6.49   0.000    -.1684304   -.0903057
gls_jstat_reti |   .4621052   .0435521    10.61   0.000     .3767446    .5474658
    gls_mstat2 |  -.0492795   .0145358    -3.39   0.001    -.0777692   -.0207899
    gls_mstat3 |  -.3209134   .0229153   -14.00   0.000    -.3658265   -.2760003
    gls_mstat4 |  -.4347647   .0213013   -20.41   0.000    -.4765145   -.3930149
    gls_mstat5 |  -.5639636   .0268586   -21.00   0.000    -.6166056   -.5113217
    gls_mstat6 |  -.2582924   .0173777   -14.86   0.000    -.2923522   -.2242327
 gls_attr_in_1 |  -.0909237   .0140381    -6.48   0.000    -.1184379   -.0634094
 gls_attr_in_2 |  -.0348121   .0130611    -2.67   0.008    -.0604114   -.0092129
 gls_attr_in_3 |  -.0768448   .0119836    -6.41   0.000    -.1003321   -.0533574
   gls_period2 |   .0014948   .0135129     0.11   0.912    -.0249901    .0279797
   gls_period3 |   .0579759   .0137876     4.20   0.000     .0309526    .0849991
   gls_period4 |  -.0510945   .0130997    -3.90   0.000    -.0767694   -.0254195
   gls_period5 |  -.0884757   .0131771    -6.71   0.000    -.1143024    -.062649
   gls_period6 |  -.0435257   .0132806    -3.28   0.001    -.0695553   -.0174961
   gls_period7 |  -.0345467   .0134656    -2.57   0.010    -.0609388   -.0081547
   gls_period8 |  -.0978498   .0137051    -7.14   0.000    -.1247112   -.0709883
   gls_period9 |  -.1856389   .0138363   -13.42   0.000    -.2127575   -.1585203
  gls_period10 |  -.1493013    .014028   -10.64   0.000    -.1767957    -.121807
       gls_one |   8.160081   1.484557     5.50   0.000     5.250402    11.06976
--------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
_all: Identity               |
   sd(gls_s~e2..gls_s~16)(1) |   .0001029   .0000373      .0000506    .0002092
-----------------------------+------------------------------------------------
                sd(Residual) |   .8962669   .0018025       .892741    .8998068
------------------------------------------------------------------------------
LR test vs. linear regression: chibar2(01) =   296.50 Prob >= chibar2 = 0.0000
(1) gls_spline2 gls_spline3 gls_spline4 gls_spline5 gls_spline6 gls_spline7 gls_spline8
    gls_spline9 gls_spline10 gls_spline11 gls_spline12 gls_spline13 gls_spline14
    gls_spline15 gls_spline16

. est store model1

. global N_fe = e(k_f)

. 
. /*  --------( prediction: smooth function )-------------------------------- */
. predict reff*, reffects level(_all) 

. global all_vars $exo_vars $spline_vars  

. foreach var of global all_vars{     
  2.     qui replace gls_`var' = `var'   
  3. }

. replace gls_one = 1
(0 real changes made)

. 
. local i = 1                 

. local j = `i' + 1

. while `i' < $N_spl {
  2. gen re_product`i' = reff`i'*spline`j'
  3. local i = `i' + 1
  4. local j = `i' + 1
  5. }

. egen blups_rcoef = rowtotal(re_product1-re_product$N_knots)

. drop re_product*

. adjust $gls_exo_vars, generate(yfixed)  

-------------------------------------------------------------------------------------------
     Dependent variable: gls_y     Equation: gls_y     Command: xtmixed
       Created variable: yfixed
   Variables left as is: gls_spline1, gls_spline~q, gls_spline~u
 Covariates set to mean: gls_ln_hhinc = 7.5732213, gls_ln_hhn~s = .91172582,
                         gls_female = .54600666, gls_bad_he~h = .60834897,
                         gls_educ_mid = .39292877, gls_educ_h~h = .37492722,
                         gls_jstat_~b = .03591415, gls_jstat_~l = .59553924,
                         gls_jstat_~m = .03246911, gls_jstat_~i = .00437504,
                         gls_mstat2 = .11654105, gls_mstat3 = .07631656,
                         gls_mstat4 = .05737692, gls_mstat5 = .01808242,
                         gls_mstat6 = .17871353, gls_attr_i~1 = .06010222,
                         gls_attr_i~2 = .09249854, gls_attr_i~3 = .11861131,
                         gls_period2 = .0814356, gls_period3 = .07960794,
                         gls_period4 = .11298279, gls_period5 = .1118587,
                         gls_period6 = .11544931, gls_period7 = .1121741,
                         gls_period8 = .10756453, gls_period9 = .10593097,
                         gls_period10 = .1036181, gls_one = 1
-------------------------------------------------------------------------------------------

----------------------
      All |         xb
----------+-----------
          |    1.84365
----------------------
     Key:  xb  =  Linear Prediction

. gen EBLUP = blups_rcoef + yfixed        

. 
. /*  --------( 95% Bonferroni-corrected CIs )------------------------------- */
. sort x

. matrix accum CTC = _gls_spline1 _gls_spline1_sq _gls_spline1_cu     ///
>     $copy_gls_exo_vars _gls_spline2-_gls_spline$N_spl, noconstant       
(obs=123656)

. global N_fe_spl = $N_fe + $N_knots

. matrix nullcol = J($N_knots,$N_fe,0)

. matrix nullrow = J($N_fe,$N_fe_spl,0)

. matrix D = (nullrow \ nullcol,I($N_knots))

. matrix estimates = e(b)                 

. matrix est_sig_e = estimates[1,"lnsig_e:_cons"]

. matrix est_sig_u = estimates[1,"lns1_1_1:_cons"]

. scalar sigma_e = exp(est_sig_e[1,1])

. scalar sigma_u = exp(est_sig_u[1,1])

. scalar lambda = (sigma_e^2/sigma_u^2) 

. matrix Cov = sigma_e^2*inv(CTC + lambda*D)

. collapse (mean) EBLUP spline1 spline1_sq spline1_cu spline2-spline$N_spl $exo_vars, by(x)

. mkmat spline1 spline1_sq spline1_cu $exo_vars spline2-spline$N_spl, matrix(C_x)

. matrix Var = C_x*Cov*C_x'

. matrix Var_f = vecdiag(Var)

. matrix Var_f = Var_f'

. svmat Var_f, names(Var_f)

. gen se_f = sqrt(Var_f)

. scalar z = invnormal(1-0.025/83)

. gen ci_upper = EBLUP + z * se_f

. gen ci_lower = EBLUP - z * se_f

. collapse (mean) EBLUP ci_*, by(x)

. 
. /*  --------( plot of smooth function )------------------------------------ */
. cap svmat knots, names(knots)

. levelsof knots, local(myknots)
21 25 29 32 35 38 41 44 48 52 56 60 65 70 77

. twoway (rarea ci_upper ci_lower x, bcolor(gs12)) ||     ///
>      (line EBLUP x, sort connect(L) lstyle(solid)) ,    ///
>      yscale(range(3.0(0.5)6)) ylabel(3.0(0.5)6)         /// 
>      legend(off) ytitle(life satisfaction) xtitle(age)  ///
>      xmtick(`myknots', tpos(inside) tlength(*2) tlc(black)) 

. log close    
